clc;
clear all;
close all;

tic % STOPWATCH ON

chann_flag = 1;
fig_on = 0;
no_noise = 0;
iter_max = 4;
inter_p = 73;
fid = fopen('result.txt', 'a');
str=datestr(clock);
fprintf(fid,'\n\n*********************%s*********************start\n', str);
% =========================================================================
sps = 8;
Rb = 2400;
Tb = 1 / Rb; % bit duration
Ts = Tb*(1/sps); 
Fs = 1 / Ts;
BW = 3600;
BW_bb = BW / 2;         % 基带带宽

%% Pulse shape & Variable initialization
pulse            = 3;      % 1 -> lorentzian pulse
                           % 2 -> GMSK pulse BT = 0.3
                           % 3 -> LRC pulse
                           % 4 -> LREC pulse
pulse_length     = 1;      % 1  -> Full response
                           % >1 -> Partial response
M_ary            = 2^1;    % M_ary symbols used 2 -> Binary
h_mod            = 0.5;    % Modulation index
[g_t,q_t]     = pulse_gen(pulse, pulse_length, 1, Fs, sps);


% d_lpf = fdesign.lowpass('Fp,Fst,Ap,Ast', 0.2, 0.3, 0.1, 60);
% % d_lpf = fdesign.lowpass('N,Fp,Fst',51, 0.1875, 0.3);
% Hd_lpf = design(d_lpf);
% % fvtool(Hd_lpf);
% h_lpf = Hd_lpf.Numerator;

load 'sh_1d3_24_128.mat';
z=24;
kb=8;
mb=16;
nb=24;
mod_odr = log2(M_ary);

%瑞利衰落信道模型
fd = 1;
delay = [0 0.00209];   % 多径时延
pd = [0 0];  
ch = rayleighchan(1/Fs, fd, delay, pd);     % 生成信道模型
ch.dopplerspectrum = doppler.gaussian;
ch.storepathgains = 1;
if (fd~=0)
    ch.storehistory = 1;
end
ch.resetbeforefiltering = 1;    % 跳频需要每次复位信道
order = 4;
mp = 2;


bit_len = kb * z;
code_len = nb * z;
syb_len = code_len / mod_odr;
tblen = 50;
rz_fg = 1;
sap_len = syb_len * sps;
sta_len = bit_len;

t_seq = 0:Ts:(sap_len+sps)*Ts-Ts;
tdt_tilt = t_seq/Tb;
ca_tilt = exp(1i*pi*h_mod*(M_ary-1)*tdt_tilt);

TREL = msk_trellis_gen(g_t, Ts, sps, h_mod, ca_tilt, fig_on);

snr_start = -2.5;
snr_step = 1;
snr_end = 20;

 
for snr = (snr_start:snr_step:snr_end)
    
    % 求解噪声功率
    eb_n0 = 10 ^ (snr/10.0);
    n0 = 1/((eb_n0)+1);
    sigma = n0 /2;
    
    be = 0;       %总的错误比特数
    fe = 0;       %总的错误帧数
    de_be = 0;    %累加可检测错误比特数
    un_be = 0;    %累加不可检测错误比特数
    de_fe = 0;    %累加不可检测错误帧数
    un_fe = 0;
    n_iter = 0;
    
    iter_flag=1;
    iter_num=0;
    
    frm_num = 0;%数据帧数
    gogo = 1;
    while gogo == 1
        bit_tx = randi([0 1], 1, bit_len);
        
        cod_tx = step(hl_enc, bit_tx.').';
        
        if (syb_len==5)
            cod_tx = [0 1 0 0 1];
        end
        
        
        %交织
        cod_tx_inter = matintrlv (cod_tx,12,48); %对信号序列的交织
%         cod_tx_inter = inter_prm(cod_tx, inter_p);
%         cod_tx_inter = cod_tx;

        
        syb_tx_ori = 1 - 2*cod_tx_inter;
        
        pos_n = length(find(syb_tx_ori==1));
        if (mod(pos_n,2)==0)
            zf = -1;
            syb_tx = [syb_tx_ori, -1];
        else
            zf = 1;
            syb_tx = [syb_tx_ori, 1];
        end



        [dat_tx, phi_tx] = msk_mod(syb_tx, g_t, Ts, sps, h_mod);
        sig_pow = var(dat_tx);
        if (fig_on)
            freq_veiw( dat_tx, 1/Ts );
        end
        
        [ dat_rx, dat_rx_e, h_0 ] = chan_rf( dat_tx, ch, fd, snr, chann_flag, order, mp, sps, sig_pow, BW_bb);
        noise_bl = dat_rx - dat_rx_e;
        snr_real = 10*log10(sig_pow/var(noise_bl));
        
%         if (no_noise)
%             dat_rx = dat_tx;
%         else
%             dat_rx = awgn(dat_tx, snr-10*log10(sps/mod_odr)+0*10*log10(nb/kb), 'measured');
% %             dat_rx = awgn(dat_tx, snr-10*log10(sps/mod_odr)+0*10*log10(nb/kb));
%         end
        
        dat_rx_tilt = dat_rx .* ca_tilt;
        
        llr_d = zeros(1, length(syb_tx));
        % 此处开始迭代     
        err_l = zeros(1, iter_max+1);
        while (iter_num<=iter_max && iter_flag~=0)

        
%             syb_rx = msk_dem_viterbi(dat_rx_tilt, TREL, tblen, rz_fg);
            [llr_rx, syb_rx] = msk_dem_logmap(dat_rx_tilt, TREL, llr_d, rz_fg);
            llr = (1/sqrt(2)) * llr_rx(1:code_len) / (8*n0);
        
%             be_l = length(find(cod_tx~=syb_rx(1:syb_len)));
%             be = be + be_l;
%             if (be_l~=0)
%                 fe = fe + 1;
%             end
%             sta_len = syb_len;
        
            % 解交织
            llr_deint = matdeintrlv (llr,12,48);     %对交织后的数据按照交织过程解交织
%             llr_deint = deint_prm(llr, inter_p);
%             llr_deint = llr;

        
            % 译码
            [ de_be, un_be, de_fe, un_fe, n_iter, cod_rx, llr_dp, c_flag, be_now] = decode_mat(...
                hl_dec, llr_deint, bit_tx, mb, nb, z, de_be, un_be, de_fe, un_fe, n_iter, 0);
%             be = de_be + un_be;
%             fe = de_fe + un_fe;
            bit_rx = cod_rx(1:bit_len);
            bit_rx_pc = (1-sign(syb_rx)) / 2;
            llr_dp = llr_dp - llr_deint;
%             llr_dp = llr_dp*0.8^iter_num;
            % 判断译码正确否，如果不正确，对llr_dp做交织，交织后重新做解调
        
            err_l(iter_num+1) = length(find(bit_rx~=bit_tx));

%             be = be + err_l;
            if (err_l(iter_num+1)~=0)
                iter_flag=1;
                llr_dd = matintrlv (llr_dp,12,48);    %交织
%                 llr_dd = inter_prm(llr_dp, inter_p);
%                 llr_dd = llr_dp;

%                 pos_nn = length(find(cod_rx==0));
%                 if (mod(pos_nn,2)==0)
%                     llr_d = [llr_dp, -100];
%                 else
%                     llr_d = [llr_dp, 100];
%                 end
                llr_d = [llr_dd, zf*100];
            else
                iter_flag=0;
            end
%             iter_flag=0;
%             iter_num=0;

            iter_num = iter_num + 1;

        end
        
        be = be + err_l(iter_num);
        if(err_l(iter_num)~=0)
            fe=fe+1;
        end
        iter_num=0;
        iter_flag=1;
        frm_num = frm_num + 1;
          
        if (mod(frm_num,10)==0)
            bit_total = sta_len * frm_num;
            ber = be / bit_total;
            fer = fe / frm_num;
            fprintf('snr=%.1f frm_num=%d ber=%.8f fer=%.8f be=%d fe=%d\n', snr, frm_num, ber, fer, be, fe);
        end
    
        if (fe>=100 && frm_num>=400)
            gogo=0;
        end
    end
    bit_total = sta_len * frm_num;
    ber = be / bit_total;
    fer = fe / frm_num;
    fprintf('snr=%.1f frm_num=%d ber=%.8f fer=%.8f be=%d fe=%d\n', snr, frm_num, ber, fer, be, fe);
    fprintf(fid, 'snr=%.1f frm_num=%d ber=%.8f fer=%.8f be=%d fe=%d\n', snr, frm_num, ber, fer, be, fe);
end

total_time = toc; % STOPWATCH OFF

